---
output:
  html_document: default
  pdf_document: default
---
# Cost-effectiveness analysis {#analysis}

```{r setup3, include=FALSE, echo = FALSE, message = FALSE, warning=FALSE}
knitr::opts_chunk$set(message = FALSE)
rm(list = ls())  # delete everything that is in R's memory
options(scipen=1000) # removal of scientific notation

path <- "C:/Users/PouwelsXGLV/Documents/Maastricht/Eribulin CEA/Model final_Gamma/" #set location of model, to adapt!

setwd(paste(path,"report", sep = "")) # set working directory to File Source

#load package
library(rms)
library(survival)
library(muhaz)
library(survminer)
library(ggplot2)
library(data.table)
library(printr)
library(dplyr)
library(flexsurv)
library(tidyverse)
library(scales)
library(matrixStats)
library(reshape2)
library(kableExtra)
library(knitr)

# load functions
source(paste(path,"R/01_fun_analyses.R", sep = ""))
source(paste(path,"R/03_fun_dowsa.R", sep = ""))
source(paste(path,"R/03_fun_scenarios.R", sep = ""))
source(paste(path,"R/03_fun_psa_inputs.R", sep = ""))
source(paste(path,"R/03_fun_threshold.R", sep = ""))

# load model inputs
load(paste(path,"data/01_model_inputs.R", sep = ""))
```

This Section presents the [deterministic results](#det-res), the [one-way determinististic sensitivity analysis results](#dowsa), the [scenario analyses results](#scenarios), the [probabilistic sensitivity analysis results](#psa) and the [threshold analysis results](#threshold) that are obtained through the cost-effectiveness model. These results are obtained by applying the functions `CEModel_det` and `CEModel_prob` that are presented in Section \@ref(model). This Section is concluded with the [value of information and risk analyses](#voi). In all these analyses, the `wtp` object, representing the willingness-to-pay threshold, has been set to &euro;80,000 per QALY.

## Deterministic discounted results {#det-res}
As explained in the [previous Section](#model), deterministic results are obtained by using the `CEModel_det` function and the `params_det` object as `params` argument of the function. The aggregated deterministic results are presented in Table \@ref(tab:res-det-disc).  
```{r, echo = T}
l.det.res <- CEModel_det(params = params_det) # run the function
```

Table: (\#tab:res-det-disc) Discounted deterministic results  

|Totals eribulin||Totals non-eribulin|| Incrementals|| ICER |  
|:------------|:-|-----------------|:-|:-------|:---|:-----|  
|QALY   |  Costs|QALY    |     Costs| QALY   |  Costs|      |    
|`r round(l.det.res[[1]][,"TQ_disc.trt"],2)`| &euro;`r format(round(l.det.res[[1]][,"TC_disc.trt"],0), big.mark = ",", scientific = FALSE)`|`r round(l.det.res[[1]][,"TQ_disc.comp"],2)`| &euro;`r format(round(l.det.res[[1]][,"TC_disc.comp"],0), big.mark = ",", scientific = FALSE)`|`r round(l.det.res[[1]][,"Inc.QALY"],2)`|&euro;`r format(round(l.det.res[[1]][,"Inc.Costs"],0), big.mark = ",", scientific = FALSE)`|&euro;`r format(round(l.det.res[[1]][,"ICER"],0), big.mark = ",", scientific = FALSE)`|

## Disaggregated deterministic undiscounted results
The deterministic undiscounted results can also be obtained from the `l.det.res` list and are provided in  as follows. Table \@ref(tab:res-det-undisc) shows the undiscounted deterministic results.  

```{r res-det-undisc, echo = F}
# we will reorganise the undiscounted output so that it looks nice
df.und.det.res <- l.det.res[[2]] # export results from the list to a separate dataframe

df.und.det.res[grep(pattern = "TC", x = rownames(df.und.det.res)),] <- paste('\u20ac',format(round(df.und.det.res[grep(pattern = "TC", x = rownames(df.und.det.res)),],0), big.mark = ",", scientific = FALSE)) # round costs to 0 decimal

df.und.det.res[grep(pattern = "TC", x = rownames(df.und.det.res), invert = TRUE),] <- round(as.numeric(as.character(df.und.det.res[grep(pattern = "TC", x = rownames(df.und.det.res), invert = TRUE),])),3) # convert LY and QALY to numbers & round LY and QALY to 3 decimals

rownames(df.und.det.res) <- gsub("TLY.", "Total life years - ", rownames(df.und.det.res))
rownames(df.und.det.res) <- gsub("TQ.", "Total QALY - ", rownames(df.und.det.res))
rownames(df.und.det.res) <- gsub("TC.", "Total costs - ", rownames(df.und.det.res))
rownames(df.und.det.res) <- gsub("pfs", " progression-free", rownames(df.und.det.res))
rownames(df.und.det.res) <- gsub("pd", " progressed disease", rownames(df.und.det.res))
rownames(df.und.det.res) <- gsub(".disutil", " disutility AE", rownames(df.und.det.res))
rownames(df.und.det.res) <- gsub(".systhx", " systemic treatment costs", rownames(df.und.det.res))
rownames(df.und.det.res) <- gsub(".adm", " administration", rownames(df.und.det.res))
rownames(df.und.det.res) <- gsub(".mon", " monitoring", rownames(df.und.det.res))
rownames(df.und.det.res) <- gsub(".surg", " surgery", rownames(df.und.det.res))
rownames(df.und.det.res) <- gsub(".hosp", " hospitalisation", rownames(df.und.det.res))
rownames(df.und.det.res) <- gsub(".radio", " radiotherapy", rownames(df.und.det.res))
rownames(df.und.det.res) <- gsub(".unrel.comp", " unrealted health care", rownames(df.und.det.res))
rownames(df.und.det.res) <- gsub(".travel", " travel", rownames(df.und.det.res))
rownames(df.und.det.res) <- gsub(".lastyear", " last year of life", rownames(df.und.det.res))
rownames(df.und.det.res) <- gsub("AE", "adverse event", rownames(df.und.det.res))
knitr::kable(
  df.und.det.res, booktabs = TRUE,
  caption = 'Undiscounted deterministic results'
)
```

## Deterministic one-way sensitivity analyses {#dowsa}
To generate the deterministic one-way sensitivity analyses (DOWSA) results, all parameters of the model were individually varied within their 95% confidence interval (CI), or by a 25% difference from their mean value when the 95%CI was not available. All calculations concerning these analyses are detailed in the  `generate_dowsa` function (`R`folder). The results of these analyses are summarised in the [tornado diagram](#fig:tornardo) below. The tornado diagram displays the 15 most influential parameters on the results. Interested readers can load the `03_dowsa_inputs.R` (`R` folder) to inspect the inputs used for the DOWSA.  

```{r tornardo, echo = TRUE, warning = FALSE, fig.cap = "Tornado diagram"}
#Draw tornado diagram
##SOURCE tornado diagram: https://stackoverflow.com/questions/55751978/tornado-both-sided-horizontal-bar-plot-in-r-with-chart-axes-crosses-at-a-given
df <- generate_dowsa(params_det = params_det)

df <- as.data.frame(df)
df$UL_Difference <- as.numeric(as.character(df$UL_Difference))
df$Lower_Bound <- as.numeric(as.character(df$Lower_Bound))
df$Upper_Bound <- as.numeric(as.character(df$Upper_Bound))
df <- df[order(df$UL_Difference, decreasing = TRUE),] #order
df <- head(df, 15) # select 15 most influential parameters

# original value of output
base.value <- as.numeric(as.character(l.det.res[[1]][1,"Inc.QALY"])) * 
  wtp - as.numeric(as.character(l.det.res[[1]][1,"Inc.Costs"])) 

# get order of parameters according to size of intervals
# (I use this to define the ordering of the factors
# which I then use to define the positions in the plot)
order.parameters <- df %>% arrange(UL_Difference) %>%
  mutate(Parameter=factor(x=Parameter, levels=Parameter)) %>%
  select(Parameter) %>% unlist() %>% levels()

# width of columns in plot (value between 0 and 1)
width <- 0.95

# get data frame in shape for ggplot and geom_rect
df.2 <- df %>% 
  # gather columns Lower_Bound and Upper_Bound into a single column using gather
  gather(key='type', value='output.value', Lower_Bound:Upper_Bound) %>%
  # just reordering columns
  select(Parameter, type, output.value, UL_Difference) %>%
  # create the columns for geom_rect
  mutate(Parameter=factor(Parameter, levels=order.parameters),
         ymin=pmin(output.value, base.value),
         ymax=pmax(output.value, base.value),
         xmin=as.numeric(Parameter)-width/2,
         xmax=as.numeric(Parameter)+width/2)
ggplot() + 
  geom_rect(data = df.2, 
            aes(ymax=ymax, ymin=ymin, xmax=xmax, xmin=xmin, fill=type)) +
  theme_bw() + 
  labs(y = "Incremental Net Monetary Benefit") +
  scale_y_continuous(labels = dollar_format(prefix = "\u20ac ", suffix = "")) +
  theme(axis.title.y=element_text(colour = "black"), legend.position = 'bottom',
        legend.title = element_blank(),
        axis.title.x = element_text(size=8)) + 
  geom_hline(yintercept = base.value) +
  scale_x_continuous(breaks = c(1:length(order.parameters)), 
                     labels = order.parameters) +
  coord_flip()
```

## Scenario analyses {#scenarios}
The following scenarios were performed during the current cost-effectiveness analysis All the calculations underlying these analyses are provided in the *03_fun_scenarios.R* file (`R` folder) and the scenario analyses results are summarised [here](#(tab:scen-res)).    
- Hospital perspective: in this scenario, all costs which were not collected within SONABRE, expect administration and preparation costs, were set to 0.  
- In the three following scenarios, we changed the distribution of OS, PFS, and TTD to Weibull, log-logistic, and lognormal respectively. This was done by modifying the `params_det$s.mod.OS`, `params_det$s.mod.PFS`, and `params_det$s.mod.TTD` objects to "Weibull", "Loglogistic", and "Lognormal" respectively.  
- Vial sharing for eribulin: in this scenario, we used the mean administered dose of eribulin (1.93) to calculate the acquisition costs of eribulin. This assumed vial sharing. This analysis was done by setting the `params_det$c.PFS.systhx.trt` to `r round(c.eri_mean,0)`, using the `c.eri_mean` object.  
- Alternative utility values from previous publication: in this scenario, the health state utility value associated with PFS (`params_det$u.PFS` object) and PD (`params_det$u.PD` object) were 0.756 and 0.496, respectively. These utility values were obtained from the evidence review group (ERG) report concerning eribulin as third (or later) treatment line in the United Kingdom. The ERG used these utility values because they assumed that the utility values reported in Lopes et al.[@Lopes2013] and used in the company submission were not correct. These corrected utility values were considered to be more approriate than the ones reported in Lopes et al [@Lopes2013].  
- No unrealted health care costs in the last year of life: in this scenario, the unrelated health care costs of the last year of life (obtained from PAID) were not accounted for. The inputs of the `params_det$c.related` object were set to 0.  
- Use of hazard ratio from EMBRACE for OS and PFS: in this scenario, the hazard ratio's reported in EMBRACE were applied on the non-eribulin survival curves to estimate PFS and OS of the eribulin group. The hazard ratios for OS and PFS were 0.81 and 0.87 respectively. This was performed by modifying the `params_det$s.OS.trt` and `params_det$s.PFS.trt` objects of the `params_det` list.  

```{r scen-res, echo = T}
df.scen <- as.data.frame(generate_scenario_results(params_det = params_det))

names(df.scen) <- c("Total QALY eribulin", "Total QALY non-eribulin", "Total costs eribulin", "Total costs non-eribulin", "Incremental Costs","Incremental QALY","ICER")
rownames(df.scen) <- c("Deterministic base-case", "Hospital perspective", 
                       "Use of Weibull distribution to model OS",
                       "Use of Weibull distribution to model PFS",
                       "Use of lognormal distribution to model TTD",
                       "Vial sharing for eribulin",
                       "Alternative utility values from ERG",
                       "No last year of life treatment",
                       "Use of hazard ratio from EMBRACE for OS and PFS")
df.scen[,c("Total costs eribulin", "Total costs non-eribulin", "Incremental Costs","ICER")] <- round(df.scen[,c("Total costs eribulin", "Total costs non-eribulin", "Incremental Costs","ICER")],0)
df.scen[,c("Total QALY eribulin", "Total QALY non-eribulin")] <- round(df.scen[,c("Total QALY eribulin", "Total QALY non-eribulin")],3)
df.scen <- df.scen[,c("Total QALY eribulin", "Total costs eribulin", "Total QALY non-eribulin", "Total costs non-eribulin", "Incremental QALY","Incremental Costs","ICER")]

knitr::kable(
  df.scen, booktabs = TRUE,
  caption = 'Scenario analyses results')
```

## Probabilistic results {#psa}
Probabilistic results are obtained by running the `CEModel_prob` function. To execute this function, probabilistic parameters have first to be estimated. Probabilistic parameters are estimated outside the "for loop" that is used within the `CEModel_prob` function. Probabilistic parameters are estimated through the `generate_psa_inputs` function with the two arguments `n.sim` (number of iterations) and `seed.num` (random seed number). In the *01_model_inputs.R* file, `n.sim` = 10,000 iterarions and `seed.num` = 500. The `generate_psa_inputs` function can be inspected in the *03_fun_psa_inputs.R* file (`R` folder). This function estimates the probabilistic parameters and combines them in a list, in this case the `l.params_prob` list. Table 2 of the manuscript provides an overview of the input parameters, their standard errors and the distribution used to estimate their probabilistic values.

```{r, echo = T, message = FALSE, warning=FALSE}
l.params_prob <- generate_psa_inputs(n.sim = n.sim, seed.num = seed.num)
``` 

Once probabilistic parameters have been estimated, one can run the `CEModel_prob` function using the code below by specifying the `params` and `n.sim` arguments. For our base-case analysis, these are respectively `l.params_prob` and `n.sim`. The results obtained from this function is a matrix, that we called `m.prob.res`. Table \@ref(tab:res-prob) shows the mean probabilistic results. Finally, the probabilistic results can be visualised in the [cost-effectiveness plane and the cost-effectiveness acceptability curve](#fig:graphs-res). The cost-effectiveness plane provides the incremental costs and QALY of each iteration while the cost-effectiveness acceptability curve provides the probability of eribulin being cost effective at different willingness-to-pay thresholds. In the cost-effectiveness plane, the dotted line represents the &euro;80,000 per QALY willingness-to-pay threshold. The functions used in this part of the code have been adapted from the `dampack` package [@R-dampack].

```{r, echo = T}
m.prob.res <- CEModel_prob(params = l.params_prob, n.sim = n.sim) # provides a matrix with the results of each iteration (10,000 in our analysis)
write.csv(m.prob.res[,1:abs(ncol(m.prob.res)-1)], file = paste(path,"output/03_probabilistic_results.csv", sep =""))# write probabilistic results in the form a of .csv file
dfprob <- as.data.frame(m.prob.res)#transform output in dataframe
```

Table: (\#tab:res-prob) Mean probabilistic results  

|Totals eribulin| |Totals non-eribulin| | Incrementals| | ICER |  
|:------------|:-|-----------------|:-|:-------|:---|:-----|  
|QALY   |  Costs|QALY    |     Costs| QALY   |  Costs|      |    
|`r round(mean(dfprob[,"TQ_disc.trt"]),3)`| &euro;`r format(round(mean(dfprob[,"TC_disc.trt"]),0), big.mark = ",", scientific = FALSE)`|`r round(mean(dfprob[,"TQ_disc.comp"]),3)`| &euro;`r format(round(mean(dfprob[,"TC_disc.comp"]),0), big.mark = ",", scientific = FALSE)`|`r round(mean(dfprob[,"Inc.QALY"]),3)`|&euro;`r format(round(mean(dfprob[,"Inc.Costs"]),0), big.mark = ",", scientific = FALSE)`|&euro;`r format(round(mean(dfprob[,"Inc.Costs"])/mean(dfprob[,"Inc.QALY"]),0), big.mark = ",", scientific = FALSE)`|  

```{r graphs-res, echo = T, fig.cap = "Cost-effectiveness plane and cost-effectiveness acceptability curve"}
limQALY  <- max(abs(dfprob$Inc.QALY)) #define limit Incremental QALY axis
limCosts <- max(abs(dfprob$Inc.Costs))#define limit Incremental Costs axis

#Cost-effectiveness plane
ggplot(dfprob, aes(x = Inc.QALY, y = Inc.Costs)) + 
  ggtitle("Cost-effectiveness plane") +
  geom_point(colour = "lightblue", shape = 1) + 
  geom_point(aes(x = mean(Inc.QALY), y = mean(Inc.Costs)), colour = 'black', shape = 18) +
  xlab ("Incremental QALY") + 
  ylab("Incremental costs") +
  geom_hline(yintercept = 0, linetype="dashed", 
             color = "red") +
  geom_vline(xintercept = 0, linetype="dashed", 
             color = "red") + 
  geom_abline(intercept = 0, slope = 80000, linetype= "dotted", 
              color = "black") + # 80,000 per QALY threshold line
  stat_ellipse(type = "t",#assumes multivariate t distribution, using "norm would assume a multivariate normal distribution, 95%CI
               colour = "blue",
               level = 0.95, size = 1.5) +
  xlim(c(-limQALY,limQALY)) +
  scale_y_continuous(labels = dollar_format(prefix = "\u20ac ", suffix = ""),
                     limits = c(-limCosts,limCosts)) +
  theme_bw()

#### Plot CEAC ####
r.CEAC <- calcCEAC(Q.trt = dfprob$TQ_disc.trt, C.trt = dfprob$TC_disc.trt, 
                   Q.comp = dfprob$TQ_disc.comp, C.comp = dfprob$TC_disc.comp, v.wtp = v.wtp)
dfCEAC <- as.data.frame(r.CEAC)
ggplot(data = dfCEAC[which(dfCEAC$WTP.threshold<110000),], aes (x= WTP.threshold, y = Prob.trt)) + 
  ggtitle("Cost-effectiveness acceptability curve") +
  geom_line(colour = "red", linejoin = "bevel", size = 1.5) + ylim(c(0, 1)) +
  xlab("Willingness to pay thresholds") + 
  scale_x_continuous(labels = dollar_format(prefix = "\u20ac ", suffix = "")) +
  ylab("Probability of being cost effective") +
  theme_bw()
```

### Threshold analysis result {#threshold}
A threshold analysis was performed to determine at which price eribulin would become cost effective. This threshold analysis was performed based on the probabilistic results and is performed through the `perform_threshold_analysis` function (*03.fun_threshold.R*, `R` folder). The argument of this function `dfprob` should be a dataframe containing the probabilistic results. The `scal.r.thr.disc` object is the threshold at which eribulin would become cost effective (the discount at which the ICER would be below the &euro;80,000 willingness-to-pay threshold).

```{r thesh, echo = T}
df.r.thresh <- perform_threshold_analysis(dfprob = dfprob)

scal.r.thr.disc <- head(df.r.thresh[which(df.r.thresh[,"ICER"] < wtp),"Discount"],1)

```

Based on this analysis, a `r scal.r.thr.disc` discount would be needed for eribulin to become cost effective.


## Value of information and risk analyses {#voi}

### Per patient results
We  calculated the expected value of perfect information (EVPI), also called Payer Uncertainty Burden (PUB), and the Payer Strategy Burden (PSB) associated with the reimbursement of eribulin [@Grimm2017], which is referred to the monetary loss associated with reimbursing eribulin the manuscript. The EVPI (or PUB) can be added to the PSB, which forms the risk per patient associated with reimbursing eribulin. To calculated the EVPI, we have adapted the `cal_evpi` function (called `evpi`) from the `dampack` package to calculate the EVPI for a range of willingness-to-pay thresholds [@R-dampack]. This function is available in the *01_fun_analysis.R* file. The EVPI per patient is plotted in Figure \@ref(fig:VOI). The risk associated with reimbursing eribulin has been estimated by calculating the difference in expected [net benefit](#glossary) between non-eribulin and eribulin (since eribulin is not cost effective). The calculations are provided below.

```{r, echo =T}
print.function(evpi)
```

``` {r VOI, echo = T, fig.cap = "Expected value of perfect information per patient"}
#EVPI
df_evpi <- evpi(v.wtp = v.wtp, m.e = cbind(dfprob$TQ_disc.trt, dfprob$TQ_disc.comp), m.c = cbind(dfprob$TC_disc.trt, dfprob$TC_disc.comp))

ggplot(data = df_evpi[which(df_evpi$WTP<110000),], aes(x = WTP , y = EVPI)) +
  #geom_point() +
  geom_line() +  
  scale_x_continuous(labels = dollar_format(prefix = "\u20ac ", suffix = ""), breaks = number_ticks(10), limits = c(0,100000)) + 
  scale_y_continuous(labels = dollar_format(prefix = "\u20ac ", suffix = ""), breaks = number_ticks(6)) +
  xlab("Willingness-to-pay (\u20ac/QALY)") +
  ylab("EVPI (\u20ac/ per patient)") +
  theme_bw(base_size = 10) +
  theme(legend.position="bottom")
```

```{r PSBPUB, echo = TRUE}
EVPI80000 <- round(df_evpi[which(df_evpi$WTP == 80000),"EVPI"],0)
PSB80000  <- round((mean(dfprob$TQ_disc.comp) * wtp - mean(dfprob$TC_disc.comp)) - (mean(dfprob$TQ_disc.trt) * wtp - mean(dfprob$TC_disc.trt)),0) #PSB at 80,000

risk80000 <- EVPI80000 + PSB80000
```
In short, the EVPI (or PUB) per patient was &euro;`r EVPI80000` and the PSB was &euro;`r PSB80000`. This resulted in a total risk  of &euro;`r risk80000` per patient. These figures correspond to `r round(EVPI80000/wtp,2)`, `r round(PSB80000/wtp,2)`, and
`r round(risk80000/wtp,2)` forgone QALY per patients.

### Population results
The above mentioned figures were scaled up to the Dutch population, assuming a technology lifetime of 10 years. The calculations for the 10-year population can be seen in the following code chunk.

```{r, echo = T}
#first calculation of number of eligible patients
scal.inci.pt   <- 14882 	#number of incident patients in 2018: https://www.iknl.nl/nkr-cijfers?fs%7Cepidemiologie_id=7&fs%7Ctumor_id=369&fs%7Cregio_id=135&fs%7Cgedrag_id=12&fs%7Cperiode_id=86%2C87%2C88%2C89%2C90%2C91%2C92%2C93%2C94%2C95%2C96%2C97%2C98%2C99%2C100%2C101%2C102%2C103%2C104%2C105%2C106%2C107%2C108%2C109%2C110%2C111%2C112%2C113%2C114%2C115&fs%7Cgeslacht_id=16&fs%7Cleeftijdsgroep_id=76%2C40%2C41%2C42%2C43%2C44%2C45&fs%7Cjaren_na_diagnose_id=18&fs%7Ceenheid_id=2&cs%7Ctype=column&cs%7CxAxis=periode_id&cs%7Cseries=leeftijdsgroep_id&ts%7CrowDimensions=periode_id&ts%7CcolumnDimensions=leeftijdsgroep_id&lang%7Clanguage=nl

scal.pr.i.abc  <- 0.05    #proportion metastatic patients at diagnosis ("incident" advanced breast cancer (ABC)  patients): https://www.zorginstituutnederland.nl/publicaties/adviezen/2017/04/11/pakketadvies-palbociclib
scal.pr.d.abc <- 0.2    #proportion of patients who will develop ABC: https://www.zorginstituutnederland.nl/publicaties/adviezen/2017/04/11/pakketadvies-palbociclib

scal.num.pt <- ceiling(scal.inci.pt*scal.pr.i.abc + scal.inci.pt*(1-scal.pr.i.abc)*scal.pr.d.abc) # total number of ABC patients

scal.pr.eli.eri <- (45+74)/721 #proportion of eligible patients for eribulin in the entire population since the introduction of eribulin

scal.eri <- 15/32 # proportion of eligible patients who actually received eribulin after 01-01-2016 (15 is the number of patients who received eribulin, 32 is the number of eligible patients)

scal.eli.pop  <- scal.num.pt*scal.pr.eli.eri # proportion of eligible patients in the first year 
scal.EVPIpop  <- ceiling(sum(sapply(1:10, function (x) scal.eli.pop/(1+params_det$d.c)^(x-1)))) #calculate size of the patient population, discounting the size of the population each year by 4%

scal.pop.EVPI <- round(scal.EVPIpop * df_evpi[which(df_evpi$WTP == 80000), "EVPI"],0) #population EVPI at the 80000 euros per QALY threshold
pop.PSB <- round(scal.EVPIpop * PSB80000 * scal.eri,0) # population PSB at the 80000 euros per QALY threshold, assuming that not all patients received eribulin

pop.EVPI.1year <- round(scal.eli.pop * df_evpi[which(df_evpi$WTP == 80000), "EVPI"],0)# yearly EVPI

pop.PSB.1year <- round(scal.eli.pop * PSB80000 * scal.eri,0)# PSB over 1 year, assuming 47% of patients would receive eribulin
```

The population results were the following. The yearly population EVPI (or PUB) was &euro;`r format(pop.EVPI.1year, big.mark=",",scientific=FALSE)` and the PSB was &euro;`r round(pop.PSB.1year/1000000,1)` million, which resulted in a risk of &euro;`r round((pop.EVPI.1year + pop.PSB.1year)/1000000,1)` million. These figures correspond to `r round(pop.EVPI.1year/wtp,2)`, `r round(pop.PSB.1year/wtp,2)`, and
`r round((pop.EVPI.1year + pop.PSB.1year)/wtp,2)` forgone QALY per year. The EVPI (or PUB) was &euro;`r format(scal.pop.EVPI, big.mark=",",scientific=FALSE)` and the PSB was &euro;`r round(pop.PSB/1000000,1)` million over a 10-year technology lifetime. This resulted in a 10-year risk of &euro;`r round((scal.pop.EVPI + pop.PSB)/1000000,1)` million.
